Let us set some global options for all code chunks in this document.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = TRUE,
# Disable warnings printed by R code chunks
warning = TRUE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
fig_count <<- fig_count + 1
paste0("Figure ", fig_count, ": ", caption)
}
# Define the function to truncate a number to two decimal places
truncate_to_two <- function(x) {
floor(x * 100) / 100
}
m1table <- rSPDE:::m1table
m2table <- rSPDE:::m2table
m3table <- rSPDE:::m3table
m4table <- rSPDE:::m4table# install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/testing"), dep=TRUE)
# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)## Loading required package: Matrix
## This is INLA_25.05.01 built 2025-05-01 18:43:33 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - List available models/likelihoods/etc with inla.list.models()
## - Use inla.doc(<NAME>) to access documentation
## - Consider upgrading R-INLA to testing[25.05.07] or stable[24.12.11].
## Loading required package: fmesher
## This is rSPDE 2.5.1
## - See https://davidbolin.github.io/rSPDE for vignettes and manuals.
## This is MetricGraph 1.4.1
## - See https://davidbolin.github.io/MetricGraph for vignettes and manuals.
##
## Attaching package: 'MetricGraph'
## The following object is masked from 'package:stats':
##
## filter
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# For each m and beta, this function returns c_m/b_{m+1} and the roots of rb and rc
my.get.roots <- function(order, beta) {
mt <- get(paste0("m", order, "table"))
rb <- rep(0, order + 1)
rc <- rep(0, order)
if(order == 1) {
rc = approx(mt$beta, mt[[paste0("rc")]], beta)$y
} else {
rc = sapply(1:order, function(i) {
approx(mt$beta, mt[[paste0("rc.", i)]], beta)$y
})
}
rb = sapply(1:(order+1), function(i) {
approx(mt$beta, mt[[paste0("rb.", i)]], xout = beta)$y
})
factor = approx(mt$beta, mt$factor, xout = beta)$y
return(list(rb = rb, rc = rc, factor = factor))
}
# Given the roots, return the polynomial coefficients in increasing order like a+bx+cx^2+...
# So if the output is c(6, -5, 1), that means 6 - 5x + x^2
poly_from_roots <- function(roots) {
coef <- 1
for (r in roots) {
coef <- convolve(coef, c(1, -r), type = "open")
}
return(coef) # returns in increasing order like a+bx+cx^2+...
}This scheme
\[\begin{equation} U^{k+1} = C^{-1}\left(\sum_{k=1}^{m+1} a_k(LC^{-1}-p_kI)^{-1} + rI\right) CU^k \end{equation}\]
compute_partial_fraction_param <- function(factor, pr_roots, pl_roots, cte) {
pr_coef <- c(0, poly_from_roots(pr_roots)) # in decreasing order like x^n+bx^(n-1)+cx^(n-2)+...
pl_coef <- poly_from_roots(pl_roots) # in decreasing order like x^n+bx^(n-1)+cx^(n-2)+...
# factor_pr_coef <- factor * pr_coef
# pr_plus_pl_coef <- factor_pr_coef + cte * pl_coef
factor_pr_coef <- pr_coef
pr_plus_pl_coef <- factor_pr_coef + cte/factor * pl_coef
res <- gsignal::residue(factor_pr_coef, pr_plus_pl_coef)
return(list(r = res$r, p = res$p, k = res$k)) # in decreasing order like x^n+bx^(n-1)+cx^(n-2)+...
}
my.fractional.operators.frac <- function(L, # kappa^2C + G
beta,
C,
scale.factor, # kappa^2
m = 1,
time_step) {
C <- Matrix::Diagonal(dim(C)[1], rowSums(C)) # lumped
Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C)) # lumped
I <- Matrix::Diagonal(dim(C)[1])
L <- L / scale.factor # C + G/kappa^2
LCi <- L %*% Ci
roots <- my.get.roots(m, beta)
poles_rs_k <- compute_partial_fraction_param(roots$factor, roots$rc, roots$rb, time_step)
# Fill list
partial.fraction.factors <- list()
for (i in 1:(m+1)) {
partial.fraction.factors[[i]] <- (LCi - poles_rs_k$p[i] * I)/poles_rs_k$r[i]
}
partial.fraction.factors[[m+2]] <- ifelse(is.null(poles_rs_k$k), 0, poles_rs_k$k) * I
return(list(Ci = Ci, C = C, LCi = LCi, L = L, m = m, beta = beta, partial.fraction.factors = partial.fraction.factors))
}
my.solver.frac <- function(obj, v){
m <- obj$m
C <- obj$C
Ci <- obj$Ci
partial.fraction.factors <- obj$partial.fraction.factors
#v <- C %*% v
output <- partial.fraction.factors[[m+2]] %*% v
for (i in 1:(m+1)) {
output <- output + solve(partial.fraction.factors[[i]], v)
}
return(Ci %*% output)
}This scheme
\[\begin{equation} U^{k+1} = \dfrac{c_m}{b_{m+1}}\dfrac{1}{\gamma}C^{-1}\prod_{k=1}^{m+1} (LC^{-1}-R_kI)^{-1}\prod_{i=1}^m (I-r_{1i}LC^{-1})\cdot CU^k \end{equation}\]
# Given a set of complex roots, this function groups complex roots into conjugate pairs
get_conjugate_pairs <- function(roots) {
used <- rep(FALSE, length(roots))
pairs <- list()
for (i in seq_along(roots)) {
if (!used[i]) {
conj_root <- Conj(roots[i])
# Find its conjugate (within tolerance)
j <- which(!used & abs(Re(roots) - Re(conj_root)) < 1e-8 & abs(Im(roots) - Im(conj_root)) < 1e-8)
j <- j[j != i]
if (length(j) > 0) {
used[c(i, j[1])] <- TRUE
pairs[[length(pairs) + 1]] <- c(roots[i], roots[j[1]])
}
}
}
return(pairs)
}
# This function builds the polynomial c_m\b_{m+1} prod_{i=1}^m (1-r_{1i}x) + tau \prod_{j=1}^{m+1} (1-r_{2j}x)
# For example compute_sum_poly(0, rep(0,2), c(2,-3,1),1) = c(6,-7,0,1) because (1-2x)(1+3x)(1-x) = 6x^3 - 7x^2 + 0x + 1
compute_sum_poly <- function(factor, pr_roots, pl_roots, cte) {
pr_coef <- c(0, poly_from_roots(pr_roots)) # in decreasing order like x^n+bx^(n-1)+cx^(n-2)+...
pl_coef <- poly_from_roots(pl_roots) # in decreasing order like x^n+bx^(n-1)+cx^(n-2)+...
pr_plus_pl_coef <- factor * pr_coef + cte * pl_coef
return(pr_plus_pl_coef) # in decreasing order like x^n+bx^(n-1)+cx^(n-2)+...
}
compute_real_roots_and_complex_coef <- function(pr_plus_pl_coef) {
pr_plus_pl_roots <- polyroot(rev(pr_plus_pl_coef)) # polyroot expects the coefficients in increasing order
real_roots <- Re(pr_plus_pl_roots[abs(Im(pr_plus_pl_roots)) < 1e-8])
complex_roots <- pr_plus_pl_roots[abs(Im(pr_plus_pl_roots)) >= 1e-8]
complex_poly_coefs <- list()
if (length(complex_roots) > 0) {
pairs <- get_conjugate_pairs(complex_roots)
for (pair in pairs) {
coef <- rev(Re(poly_from_roots(c(pair[1], pair[2])))) # returns in x^2 + bx + c order
complex_poly_coefs[[length(complex_poly_coefs) + 1]] <- round(coef, 12)
}
} else {
complex_poly_coefs <- list() # No complex roots
}
return(list(new_factor = pr_plus_pl_coef[1],
pr_plus_pl_roots = pr_plus_pl_roots,
real_roots = real_roots,
complex_poly_coefs = complex_poly_coefs))
}
my.fractional.operators.rat <- function(L, # kappa^2C + G
beta,
C,
scale.factor, # kappa^2
m = 1,
time_step) {
C <- Matrix::Diagonal(dim(C)[1], rowSums(C)) # lumped
Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C)) # lumped
I <- Matrix::Diagonal(dim(C)[1])
L <- L / scale.factor # C + G/kappa^2
LCi <- L %*% Ci
roots <- my.get.roots(m, beta)
Pl.roots <- roots$rb
Pr.roots <- roots$rc
factor <- roots$factor
new_factor_and_roots <- compute_real_roots_and_complex_coef(compute_sum_poly(factor, Pr.roots, Pl.roots, time_step))
new_real_roots <- new_factor_and_roots$real_roots
new_factor <- new_factor_and_roots$new_factor
complex_poly_coefs <- new_factor_and_roots$complex_poly_coefs
# Fill Pr_plus_Pl.factors
Pr_plus_Pl.factors <- list()
if (length(new_real_roots) >= 1) {
for (i in 1:length(new_real_roots)) {
Pr_plus_Pl.factors[[i]] <- LCi - new_real_roots[i] * I
}
}
length_now <- length(Pr_plus_Pl.factors)
if(length(complex_poly_coefs) >= 1) {
for (i in 1:length(complex_poly_coefs)) {
}
Pr_plus_Pl.factors[[length_now + i]] <- LCi %*% LCi + complex_poly_coefs[[i]][2] * LCi + complex_poly_coefs[[i]][3] * I
}
# Fill Pr.factors
Pr.factors <- list()
Pr.factors[[1]] <- I - LCi * roots$rc[1]
if (length(roots$rc) > 1) {
for (i in 2:length(roots$rc)) {
Pr.factors[[i]] <- I - LCi * roots$rc[i]
}
}
return(list(
Ci = Ci, C = C, LCi = LCi, L = L, m = m, beta = beta,
factor = factor,
new_factor = new_factor,
Pr.factors = Pr.factors,
Pr_plus_Pl.factors = Pr_plus_Pl.factors
))
}
my.solver.rat <- function(obj, v){
Pr.factors <- obj$Pr.factors
Pr_plus_Pl.factors <- obj$Pr_plus_Pl.factors
m <- obj$m
C <- obj$C
Ci <- obj$Ci
factor <- obj$factor
new_factor <- obj$new_factor
if (m==1){
temp <- solve(Pr_plus_Pl.factors[[2]],
Pr.factors[[1]] %*% solve(
Pr_plus_Pl.factors[[1]],
C%*%v))
} else if (m==2){
temp <- solve(Pr_plus_Pl.factors[[3]],
Pr.factors[[2]] %*% solve(
Pr_plus_Pl.factors[[2]],
Pr.factors[[1]] %*% solve(
Pr_plus_Pl.factors[[1]],
C%*%v)))
} else if (m==3){
temp <- solve(Pr_plus_Pl.factors[[4]],
Pr.factors[[3]] %*% solve(
Pr_plus_Pl.factors[[3]],
Pr.factors[[2]] %*% solve(
Pr_plus_Pl.factors[[2]],
Pr.factors[[1]] %*% solve(
Pr_plus_Pl.factors[[1]],
C%*%v))))
} else if (m==4){
temp <- solve(Pr_plus_Pl.factors[[5]],
Pr.factors[[4]] %*% solve(
Pr_plus_Pl.factors[[4]],
Pr.factors[[3]] %*% solve(
Pr_plus_Pl.factors[[3]],
Pr.factors[[2]] %*% solve(
Pr_plus_Pl.factors[[2]],
Pr.factors[[1]] %*% solve(
Pr_plus_Pl.factors[[1]],
C%*%v)))))
}
return((factor/new_factor) * Ci %*% temp)
}We want to solve the fractional diffusion equation \[\begin{equation} \label{eq:maineq} \partial_t u+(\kappa^2-\Delta_\Gamma)^{\alpha/2} u=f \text { on } \Gamma \times(0, T), \quad u(0)=u_0 \text { on } \Gamma, \end{equation}\] where \(u\) satisfies the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \left\{\phi\in C(\Gamma)\;\Big|\; \forall v\in V: \sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\} \end{equation}\]
If \(f=0\), then the solution is given by \[\begin{equation} \label{eq:sol_reprentation} u(s,t) = \displaystyle\sum_{j\in\mathbb{N}}e^{-\lambda^{\alpha/2}_jt}\left(u_0, e_j\right)_{L_2(\Gamma)}e_j(s). \end{equation}\]
# Function to build a tadpole graph and create a mesh
gets_graph_tadpole <- function(h){
edge1 <- rbind(c(0,0),c(1,0))
theta <- seq(from=-pi,to=pi,length.out = 10000)
edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
edges = list(edge1, edge2)
graph <- metric_graph$new(edges = edges)
graph$set_manual_edge_lengths(edge_lengths = c(1,2))
graph$build_mesh(h = h)
return(graph)
}Let \(\Gamma_T = (\mathcal{V},\mathcal{E})\) characterize the tadpole graph with \(\mathcal{V}= \{v_1,v_2\}\) and \(\mathcal{E}= \{e_1,e_2\}\). The left edge \(e_1\) has length 1 and the circular edge \(e_2\) has length 2. As discussed before, a point on \(e_1\) is parameterized via \(s=\left(e_1, t\right)\) for \(t \in[0,1]\) and a point on \(e_2\) via \(s=\left(e_2, t\right)\) for \(t\in[0,2]\). One can verify that \(-\Delta_\Gamma\) has eigenvalues \(0,\left\{(i \pi / 2)^2\right\}_{i \in \mathbb{N}}\) and \(\left\{(i \pi / 2)^2\right\}_{2 i \in \mathbb{N}}\) with corresponding eigenfunctions \(\phi_0\), \(\left\{\phi_i\right\}_{i \in \mathbb{N}}\), and \(\left\{\psi_i\right\}_{2 i \in \mathbb{N}}\) given by \(\phi_0(s)=1 / \sqrt{3}\) and \[\begin{equation*} \phi_i(s)=C_{\phi, i}\begin{cases} -2 \sin (\frac{i\pi}{2}) \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \sin (i \pi t / 2), & s \in e_2, \end{cases}, \quad \psi_i(s)=\frac{\sqrt{3}}{\sqrt{2}} \begin{cases} (-1)^{i / 2} \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \cos (\frac{i \pi t}{2}), & s \in e_2, \end{cases}, \end{equation*}\] where \(C_{\phi, i}=1\) if \(i\) is even and \(C_{\phi, i}=1 / \sqrt{3}\) otherwise. Moreover, these functions form an orthonormal basis for \(L_2(\Gamma_T)\).
# Function to compute the eigenfunctions
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2])
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2])
if(k==0){
f.e1 <- rep(1,length(x1))
f.e2 <- rep(1,length(x2))
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2)
f.e2 <- sin(pi*k*x2/2)
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
if((k %% 2)==1){
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
f.e2 <- cos(pi*k*x2/2)
f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1])
f <- list(phi=f1,psi=f2/sqrt(3/2))
}
}
return(f)
}Let \(\alpha\in(0,2]\) and \(U_h^\tau\) denote the sequence of approximations of the solution to the weak form of problem \(\eqref{eq:maineq}\) at each time step on a mesh indexed by \(h\). Let \(z=0\) and \(U^0_h = P_hu_0\). For \(k=0,\dots, N-1\), \(U_h^{k+1}\in V_h\) solves the following scheme \[\begin{align} \label{system:fully_discrete_scheme} \langle\delta U_h^{k+1},\phi\rangle + \mathfrak{a}(U_h^{k+1},\phi) = \langle f^{k+1},\phi\rangle ,\quad\forall\phi\in V_h, \end{align}\] where \(f^{k+1} = \displaystyle\dfrac{1}{\tau}\int_{t_k}^{t^{k+1}}f(t)dt\). The solution can be represented as \[\begin{align*} U_h^k(s) = \sum_{j=1}^{N_h}u_j^k\psi_j(s) \end{align*}\] Replacing this into \(\eqref{system:fully_discrete_scheme}\) yields the following linear system \[\begin{align*} \sum_{j=1}^{N_h}u_j^{k+1}[(\psi_j,\psi_i)_{L_2(\Gamma)}+ \tau\mathfrak{a}(\psi_j,\psi_i)] = \sum_{j=1}^{N_h}u_j^{k}(\psi_j,\psi_i)_{L_2(\Gamma)}+\tau( f^{k+1},\psi_i)_{L_2(\Gamma)} \end{align*}\] for \(i = 1,\dots, N_h\). In matrix notation, \[\begin{align} \label{diff_eq_discrete} (C+\tau L^{\alpha/2})U^{k+1} = CU^k+\tau F^{k+1}, \end{align}\] where \(C\) has entries \(C_{ij} = (\psi_j,\psi_i)_{L_2(\Gamma)}\), \(L^{\alpha/2}\) has entries \(\mathfrak{a}(\psi_j,\psi_i)\), \(U^k\) has entries \(u_j^k\), and \(F^k\) has entries \(( f^{k},\psi_i)_{L_2(\Gamma)}\). For simplicity, let \(f=0\) and so we arrive at
\[\begin{equation} (C+\tau L^{\alpha/2})U^{k+1} = CU^{k} \end{equation}\]
Apply \(L^{-\alpha/2}\) to both sides to get
\[\begin{equation} (L^{-\alpha/2}C+\tau I)U^{k+1} = L^{-\alpha/2}CU^{k} \end{equation}\]
Approximate \(L^{-\alpha/2}\) by \(P_rP_l^{-1}\) to arrive at
\[\begin{equation} (P_rP_l^{-1}C+\tau I)U^{k+1} = P_rP_l^{-1}CU^{k} \end{equation}\]
Actually, \(L^{-\alpha/2}\) should be approximated by \(P_l^{-T}P_r^T\) (this is like this because we want to use \(P_r\) and \(P_l\) the way they where constructed) and so we can write
\[\begin{equation} (P_l^{-T}P_r^TC+\tau I)U^{k+1} = P_l^{-T}P_r^TCU^{k} \end{equation}\]
From this, we arrive at the scheme
\[\begin{equation} (P_r^TC+\tau P_\ell^T)U^{k+1} = P_r^TCU^k \label{eq:scheme} \end{equation}\] where (here we are using \(P_r\) and \(P_l\) the way they were constructed) \[\begin{equation} P_r = c_m\prod_{i=1}^m (I-r_{1i}C^{-1}L)\quad\text{and}\quad P_\ell = b_{m+1}C\prod_{j=1}^{m+1} (I-r_{2j}C^{-1}L) \end{equation}\] Observe that \[\begin{equation} P_r^T = c_m\prod_{i=1}^m (I-r_{1i}LC^{-1})\quad\text{and}\quad P_\ell^T = b_{m+1}\prod_{j=1}^{m+1} (I-r_{2j}LC^{-1})\cdot C \end{equation}\] since \(L\) and \(C^{-1}\) are symmetric and the factors in the product commute (because the factors in the product that make up a polynomial commute, for example, \((1-ax)(1-bx) = (1-bx)(1-ax)\)). Replacing these two in our scheme we get \[\begin{equation} \left(\dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})+\tau \prod_{j=1}^{m+1} (I-r_{2j}LC^{-1})\right)CU^{k+1} = \dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})\cdot CU^k \end{equation}\]
We now need to compute the roots of the polynomial \[\begin{equation} R(x) = \dfrac{c_m}{b_{m+1}} \prod_{i=1}^m (1-r_{1i}x)+\tau \prod_{j=1}^{m+1} (1-r_{2j}x) \end{equation}\] \(R\) has leading coefficient \(\gamma=\tau (-1)^{m+1} \prod_{j=1}^{m+1}r_{2j}\) and \(m+1\) roots \(R_k\) for \(k=1,\ldots,m+1\). That is, \[\begin{equation} R(x) = \tau(-1)^{m+1} \prod_{j=1}^{m+1}r_{2j}\prod_{k=1}^{m+1} (x-R_k)= \gamma\prod_{k=1}^{m+1} (x-R_k) \end{equation}\] We can then write our scheme as \[\begin{equation} \left(\gamma\prod_{k=1}^{m+1} (LC^{-1}-R_kI)\right)CU^{k+1} = \dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})\cdot CU^k \end{equation}\] That is, \[\begin{equation} U^{k+1} = \dfrac{c_m}{b_{m+1}}\dfrac{1}{\gamma}C^{-1}\prod_{k=1}^{m+1} (LC^{-1}-R_kI)^{-1}\prod_{i=1}^m (I-r_{1i}LC^{-1})\cdot CU^k \end{equation}\]
Roots \(r_{1i}\) and \(r_{2j}\) are real but \(R_k\) can be complex. In that case, we consider the quadratic terms. For example, if \(m=3\), and \(R\) has one real root \(r\) and two complex conjugate roots \(z\) and \(\bar{z}\) such that \((x-z)(x-\bar{z}) = x^2+ax+b\), we will work with \[\begin{equation} R(x) = \gamma(x^2+ax+b)(x-r) \end{equation}\] instead of \[\begin{equation} R(x) = \gamma(x-z)(x-\bar{z})(x-r) \end{equation}\]
From the scheme
\[\begin{equation} \left(\dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})+\tau \prod_{j=1}^{m+1} (I-r_{2j}LC^{-1})\right)CU^{k+1} = \dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})\cdot CU^k \end{equation}\]
we write
\[\begin{equation} U^{k+1} = C^{-1}\left(\dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})+\tau \prod_{j=1}^{m+1} (I-r_{2j}LC^{-1})\right)^{-1}\dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (I-r_{1i}LC^{-1})\cdot CU^k \end{equation}\]
We want the partial fraction decomposition of
\[\begin{equation} \dfrac{\dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (1-r_{1i}x)}{\dfrac{c_m}{b_{m+1}}\prod_{i=1}^m (1-r_{1i}x)+\tau \prod_{j=1}^{m+1} (1-r_{2j}x)} \end{equation}\]
or
\[\begin{equation} \dfrac{\prod_{i=1}^m (1-r_{1i}x)}{\prod_{i=1}^m (1-r_{1i}x)+\tau \dfrac{b_{m+1}}{c_m} \prod_{j=1}^{m+1} (1-r_{2j}x)} \end{equation}\]
That is
\[\begin{equation} \sum_{k=1}^{m+1} \dfrac{a_k}{x-p_k} + r = \sum_{k=1}^{m+1} a_k(x-p_k)^{-1} + r \end{equation}\]
Therefore
\[\begin{equation} U^{k+1} = C^{-1}\left(\sum_{k=1}^{m+1} a_k(LC^{-1}-p_kI)^{-1} + rI\right) CU^k \end{equation}\]
Therefore
\[\begin{equation} U^{k+1} = C^{-1}\left(\sum_{k=1}^{m+1} a_k(LC^{-1}-p_kI)^{-1} + rI\right) (CU^k+\tau F^{k+1}) \end{equation}\]
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.16 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
T_final <- 2
time_step <- 0.01
time_seq <- seq(0, T_final, by = time_step)
# Compute the FEM matrices
graph$compute_fem()
G <- graph$mesh$G
C <- graph$mesh$C
I <- Matrix::Diagonal(nrow(C))
x <- graph$mesh$V[, 1]
y <- graph$mesh$V[, 2]
edge_number <- graph$mesh$VtE[, 1]
pos <- sum(edge_number == 1)+1
order_to_plot <- function(v)return(c(v[1], v[3:pos], v[2], v[(pos+1):length(v)], v[2]))
weights <- graph$mesh$weights# Parameters to construct U_0
N_finite <- 4 # choose an even number
adjusted_N_finite <- N_finite + N_finite/2 + 1
EIGENVAL_ALPHA <- NULL
EIGENFUN <- NULL
INDEX <- NULL
PHI_OR_PSI <- NULL
for (j in 0:N_finite) {
lambda_j_alpha <- (kappa^2 + (j*pi/2)^2)^(alpha/2)
e_j <- tadpole.eig(j,graph)$phi
EIGENVAL_ALPHA <- c(EIGENVAL_ALPHA, lambda_j_alpha)
EIGENFUN <- cbind(EIGENFUN, e_j)
INDEX <- c(INDEX, j)
PHI_OR_PSI <- c(PHI_OR_PSI, "phi")
if (j>0 && (j %% 2 == 0)) {
lambda_j_alpha <- (kappa^2 + (j*pi/2)^2)^(alpha/2)
e_j <- tadpole.eig(j,graph)$psi
EIGENVAL_ALPHA <- c(EIGENVAL_ALPHA, lambda_j_alpha)
EIGENFUN <- cbind(EIGENFUN, e_j)
INDEX <- c(INDEX, j)
PHI_OR_PSI <- c(PHI_OR_PSI, "psi")
}
}
# Building the initial condition as \sum coeff_j EIGENFUN_j
coeff <- 50*(1:adjusted_N_finite)^-1
coeff[-5] <- 0
U_0 <- EIGENFUN %*% coeff
# Building the true solution as \sum coeff_j EIGENFUN_j e^{-\lambda_j^{\frac{\alpha}{2}}t}
U_true <- matrix(NA, nrow = length(x), ncol = length(time_seq))
for (k in 1:length(time_seq)) {
aux_k <- rep(0, length(x))
for (j in 1:adjusted_N_finite) {
aux_k <- aux_k + exp(-time_seq[k]*EIGENVAL_ALPHA[j])*coeff[j]*EIGENFUN[, j]
}
U_true[, k] <- aux_k
}
c_k <- 10
what_eigenfunction_for_ff <- 7
ff <- function(t){
return(c_k*EIGENFUN[,what_eigenfunction_for_ff]*exp(-t*EIGENVAL_ALPHA[what_eigenfunction_for_ff]))
}
FF_true <- matrix(NA, nrow = length(x), ncol = length(time_seq))
FF_sol_true <- matrix(NA, nrow = length(x), ncol = length(time_seq))
for (k in 1:length(time_seq)) {
FF_true[, k] <- ff(time_seq[k]) # this is the right hand side function
FF_sol_true[, k] <- time_seq[k]*FF_true[, k] # this is the second term in the solution
}
U_true <- U_true + FF_sol_true
graph_finer <- gets_graph_tadpole(h = 0.001)## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.15 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
finer_loc <- graph_finer$get_mesh_locations()
A <- graph$fem_basis(finer_loc)
graph_finer$compute_fem()
finer_weights <- graph_finer$mesh$weights
finer_C <- graph_finer$mesh$C
EIGENFUN_FOR_FF <- tadpole.eig(INDEX[what_eigenfunction_for_ff], graph_finer)
if (PHI_OR_PSI[what_eigenfunction_for_ff] == "phi"){
eigenfun_for_ff <- EIGENFUN_FOR_FF$phi
} else if (PHI_OR_PSI[what_eigenfunction_for_ff] == "psi"){
eigenfun_for_ff <- EIGENFUN_FOR_FF$psi
}
int_basis_eigen <- as.vector(t(as.matrix(eigenfun_for_ff)) %*% finer_C %*% A) #sum(finer_weights*eigenfun_for_ff*A[,10])
COEF <- c_k*exp(-time_seq*EIGENVAL_ALPHA[what_eigenfunction_for_ff])
FF_approx <- int_basis_eigen %*% t(COEF){r}
# Initial condition
U_0 <- 10*exp(-((x-1)^2 + (y)^2))
U_true <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_true[, 1] <- U_0
n_finite <- 1/h
for (k in 1:(length(time_seq) - 1)) {
aux_k <- rep(0, nrow(C))
for (j in 0:n_finite) {
decay_j <- exp(-time_seq[k+1]*(kappa^2 + (j*pi/2)^2)^(alpha/2))
e_j <- tadpole.eig(j,graph)$phi
aux_k <- aux_k + decay_j*sum(U_0*e_j*weights)*e_j
if (j>0 && (j %% 2 == 0)) {
e_j <- tadpole.eig(j,graph)$psi
aux_k <- aux_k + decay_j*sum(U_0*e_j*weights)*e_j
}
}
U_true[, k + 1] <- aux_k
}
{r}
my_op_rat <- my.fractional.operators.rat(L, beta, C, scale.factor = kappa^2, m = m, time_step)
U_approx1 <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_approx1[, 1] <- U_0
# Time-stepping loop
for (k in 1:(length(time_seq) - 1)) {
U_approx1[, k + 1] <- as.matrix(my.solver.rat(my_op_rat, U_approx1[, k]))
}
my_op_frac <- my.fractional.operators.frac(L, beta, C, scale.factor = kappa^2, m = m, time_step)
U_approx2 <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_approx2[, 1] <- U_0
# Time-stepping loop
for (k in 1:(length(time_seq) - 1)) {
U_approx2[, k + 1] <- as.matrix(my.solver.frac(my_op_frac, my_op_frac$C %*% U_approx2[, k] + time_step * FF_approx[, k + 1]))
}{r}
op <- fractional.operators(L, beta, C, scale.factor = kappa^2, m = m)
Pl <- op$Pl
Pr <- op$Pr
Ci <- op$Ci
C <- op$C
LHS <- t(Pr)%*% C + time_step * t(Pl)
# Initialize U matrix to store solution at each time step
U_approx2 <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_approx2[, 1] <- U_0
# Time-stepping loop
for (k in 1:(length(time_seq) - 1)) {
# Compute the right-hand side for the second equation
RHS <- t(Pr)%*% C %*% U_approx2[, k]
U_approx2[, k + 1] <- as.matrix(solve(LHS, RHS))
}
{r}
op <- fractional.operators(L, beta, C, scale.factor = kappa^2, m = m)
Pl <- op$Pl
Pr <- op$Pr
Ci <- op$Ci
C <- op$C
Pr.apply.mult <- function(v){return(Pr.mult(op, v))}
Pl.apply.solve <- function(v){return(Pl.solve(op, v))}
PliC <- apply(C, 2, Pl.apply.solve) # PlC^-1
# Precompute the LHS1 matrix
aux <- apply(PliC, 2, Pr.apply.mult) #PrPl^-1C
LHS <- aux + time_step * Matrix::Diagonal(nrow(C))
# Initialize U matrix to store solution at each time step
U_approx2 <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_approx2[, 1] <- U_0
# Time-stepping loop
for (k in 1:(length(time_seq) - 1)) {
# Compute the right-hand side for the second equation
RHS <- aux %*% U_approx2[, k]
U_approx2[, k + 1] <- as.matrix(solve(LHS, RHS))
}
x <- order_to_plot(x)
y <- order_to_plot(y)
#U_approx2 <- U_approx1
U_approx1 <- U_approx2
mean_w <- function(v){return(mean(v*weights))}
max_error_at_each_time1 <- apply((U_true - U_approx1)^2, 2, mean)
max_error_at_each_time2 <- apply((U_true - U_approx2)^2, 2, mean)
max_error_between_both_approx <- apply((U_approx1 - U_approx2)^2, 2, mean)
error1 <- sqrt(time_step * sum(max_error_at_each_time1))
error2 <- sqrt(time_step * sum(max_error_at_each_time2))
errorb <- sqrt(time_step * sum(max_error_between_both_approx))
U_true <- apply(U_true, 2, order_to_plot)
U_approx1 <- apply(U_approx1, 2, order_to_plot)
U_approx2 <- apply(U_approx2, 2, order_to_plot)
# Create interactive plot
fig <- plot_ly()
# Add first line (max_error_at_each_time1)
fig <- fig %>% add_trace(
x = ~time_seq, y = ~max_error_at_each_time1, type = 'scatter', mode = 'lines+markers',
line = list(color = 'red', width = 2),
marker = list(color = 'red', size = 4),
name = paste0("Error True and Approx 1: ", sprintf("%.3e", error1))
)
# Add second line (max_error_at_each_time2)
fig <- fig %>% add_trace(
x = ~time_seq, y = ~max_error_at_each_time2, type = 'scatter', mode = 'lines+markers',
line = list(color = 'blue', width = 2),
marker = list(color = 'blue', size = 4),
name = paste0("Error True and Approx 2: ", sprintf("%.3e", error2))
)
# Add third line (max_error_between_both_approx)
fig <- fig %>% add_trace(
x = ~time_seq, y = ~max_error_between_both_approx, type = 'scatter', mode = 'lines+markers',
line = list(color = 'orange', width = 2),
marker = list(color = 'orange', size = 4),
name = paste0("Error Between Approximations: ", sprintf("%.3e", errorb))
)
# Layout
fig <- fig %>% layout(
title = "Error at Each Time Step",
xaxis = list(title = "Time"),
yaxis = list(title = "Error"),
legend = list(x = 0.1, y = 0.9)
)
plot_data <- data.frame(
x = rep(x, times = ncol(U_true)),
y = rep(y, times = ncol(U_true)),
z_true = as.vector(U_true),
z_approx1 = as.vector(U_approx1),
z_approx2 = as.vector(U_approx2),
frame = rep(time_seq, each = length(x))
)
# Compute axis limits
x_range <- range(x)
y_range <- range(y)
z_range <- range(c(U_true, U_approx1, U_approx2))
# Initial plot setup (first frame only)
p <- plot_ly(plot_data, frame = ~frame) %>%
add_trace(
x = ~x, y = ~y, z = ~z_true,
type = "scatter3d", mode = "lines",
name = "True",
line = list(color = "green", width = 2)
) %>%
add_trace(
x = ~x, y = ~y, z = ~z_approx1,
type = "scatter3d", mode = "lines",
name = "Approx 1",
line = list(color = "red", width = 2)
) %>%
add_trace(
x = ~x, y = ~y, z = ~z_approx2,
type = "scatter3d", mode = "lines",
name = "Approx 2",
line = list(color = "blue", width = 2)
) %>%
layout(
scene = list(
xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "Value", range = z_range),
aspectratio = list(x = 2.4, y = 1.2, z = 1.2),
camera = list(
eye = list(x = 1.5, y = 1.5, z = 1), # Adjust the viewpoint
center = list(x = 0, y = 0, z = 0))
),
updatemenus = list(
list(
type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 100, redraw = TRUE), fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
)
)
),
title = "Time: 0"
)
# Convert to plotly object with frame info
pb <- plotly_build(p)
# Inject custom titles into each frame
for (i in seq_along(pb$x$frames)) {
t <- time_seq[i]
err <- signif(max_error_between_both_approx[i], 4)
pb$x$frames[[i]]$layout <- list(title = paste0("Time: ", t, " | Error: ", err))
}## This is to plot the right hand side alone
FF_true <- apply(FF_true, 2, order_to_plot)
FF_sol_true <- apply(FF_sol_true, 2, order_to_plot)
FF_approx <- apply(FF_approx, 2, order_to_plot)
plot_data <- data.frame(
x = rep(x, times = ncol(FF_true)),
y = rep(y, times = ncol(FF_true)),
ff_true = as.vector(FF_true),
ff_sol_true = as.vector(FF_sol_true),
ff_approx = as.vector(FF_approx),
frame = rep(time_seq, each = length(x))
)
# Compute axis limits
x_range <- range(x)
y_range <- range(y)
z_range <- range(c(FF_true, FF_sol_true, FF_approx))
# Initial plot setup (first frame only)
p_ff <- plot_ly(plot_data, frame = ~frame) %>%
add_trace(
x = ~x, y = ~y, z = ~ff_true,
type = "scatter3d", mode = "lines",
name = "f(s,t)",
line = list(color = "green", width = 2)
) %>%
add_trace(
x = ~x, y = ~y, z = ~ff_sol_true,
type = "scatter3d", mode = "lines",
name = "tf(s,t) = u(s,t)-SOL(u_0)",
line = list(color = "red", width = 2)
) %>%
add_trace(
x = ~x, y = ~y, z = ~ff_approx,
type = "scatter3d", mode = "lines",
name = "F^k = (f^k, phi)",
line = list(color = "blue", width = 2)
) %>%
layout(
scene = list(
xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "Value", range = z_range),
aspectratio = list(x = 2.4, y = 1.2, z = 1.2),
camera = list(
eye = list(x = 1.5, y = 1.5, z = 1), # Adjust the viewpoint
center = list(x = 0, y = 0, z = 0))
),
updatemenus = list(
list(
type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 100, redraw = TRUE), fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
)
)
),
title = "Time: 0"
)
# Convert to plotly object with frame info
pb_ff <- plotly_build(p_ff)
# Inject custom titles into each frame
for (i in seq_along(pb_ff$x$frames)) {
t <- time_seq[i]
pb_ff$x$frames[[i]]$layout <- list(title = paste0("Time: ", t))
}